/* Simulations of reallocating panels at the intrastate level.

NOTES: All directories and paths should be set in the "Reallocation_master.do" file.

*/

clear all

import delimited `"${data_path}/Generation_and_damages_withCO2_2019-11-24_eGrid_SR.csv"', case(preserve) 
rename zip zipcode
drop num_systems

merge 1:1 zipcode using `"${data_intermediate}/NumberSystems.dta"'
*Note: We have 1478 zipcodes with solar panels but no damage calculations.
*So I assign damages to these zipcodes based on the damages from the previous zipcode (when zipcodes are listed in numerical order)
replace total_damages = total_damages[_n-1] if _merge==2
drop _merge

merge 1:1 zipcode using `"${data_path}/HousingUnits"'
*need to drop zipcodes with missing housing units or zip codes that have 0 housing units
*Note: There are 121 zip codes with a positive number of installations but 0 housing units
*This can be observed by uncommenting the next line.
*count if single_units ==0 & num_systems>0 & !missing(num_systems)
drop if missing(single_units) | single_units==0

rename _merge _merge_housing

*fill in missing state variable with ZipCodeToState file
merge 1:1 zipcode using `"${data_path}/ZipCodesToState"'
replace state = stateabbreviation if missing(state)

*label the Puerto Rico zipcodes
replace state="PR" if missing(state) & zipcode>600 & zipcode<1000

*drop Puerto Rico, Hawaii, and Alaska
drop if state=="HI" | state=="AK" | state=="PR"
*drop zips if we have no data for them
drop if _merge==2
drop _merge

*Note: We have 2177 zipcodes with single occupancy housing units but no damage calculations.
*So I assign damages to these zipcodes based on the damages from the previous zipcode (when zipcodes are listed in numerical order)
replace total_damages = total_damages[_n-1] if missing(total_damages)

drop single_units_MOE _merge_housing

*generate total damages (damages per 4kW installation * number of installations)
gen total_damages_all = total_damages*num_systems

sort total_damages single_units
egen rank=rank(-total_damages)
sort rank single_units
drop rank
gen rank=_n /*ranks zipcodes by total_damages and number of houses*/

sort state rank
bysort state: egen rankbystate=rank(rank)
bysort state (rankbystate): egen sumsystemsState = sum(num_systems)
egen stategroup = group(state)
gen id = _n /*Create new overall ranking for purpose of calling values later*/

* Simulation 1 (moving all panels to highest total_damages zipcode in state subject to #of housing units)
*drop remaining
gen state_new_num_systems=0
gen remaining = sumsystemsState

sum stategroup, meanonly
local numstates=r(max)
forvalues i=1/`numstates' {
	sum rankbystate if stategroup==`i', meanonly
	local maxRankInState=r(max)
	forvalues j = 1/`maxRankInState' {
		sum id if stategroup==`i' & rankbystate==`j', meanonly
		local index1= int(r(min))
		quietly replace state_new_num_systems=min(remaining[`index1'],single_units[`index1']) if id==`index1' & remaining[`index1']>0
		quietly drop remaining
		quietly bysort state (rankbystate): egen sumNewByState = sum(state_new_num_systems)
		quietly gen remaining=sumsystemsState-sumNewByState
		quietly drop sumNewByState
	}
}

gen newState_total_damages=state_new_num_systems*total_damages
* check
bysort state (rankbystate): egen sumNewByState = sum(state_new_num_systems)
egen sumNew=sum(state_new_num_systems)
egen SumNewDamagesState=sum(newState_total_damages)
egen SumOldDamages=sum(total_damages_all)

* Simulation 2 (moving all panels to highest total_damages zipcode in state subject to 30% x #of housing units)
drop remaining
gen single_units30 = floor(0.30*single_units)
gen state_new_num_systems30=0
gen remaining = sumsystemsState

sum stategroup, meanonly
local numstates=r(max)
forvalues i=1/`numstates' {
	sum rankbystate if stategroup==`i', meanonly
	local maxRankInState=r(max)
	forvalues j = 1/`maxRankInState' {
		sum id if stategroup==`i' & rankbystate==`j', meanonly
		local index1= int(r(min))
		quietly replace state_new_num_systems30=min(remaining[`index1'],single_units30[`index1']) if id==`index1' & remaining[`index1']>0
		quietly drop remaining
		quietly bysort state (rankbystate): egen sumNewByState30 = sum(state_new_num_systems30)
		quietly gen remaining=sumsystemsState-sumNewByState30
		quietly drop sumNewByState30
	}
}

gen newState_total_damages30=state_new_num_systems30*total_damages
* check
bysort state (rankbystate): egen sumNewByState30 = sum(state_new_num_systems30)
egen sumNew30=sum(state_new_num_systems30)
egen SumNewDamagesState30=sum(newState_total_damages30)

*Calculate percentage gains in total avoided damages
gen gain = (SumNewDamagesState - SumOldDamages)/SumOldDamages
gen gain30 = (SumNewDamagesState30 - SumOldDamages)/SumOldDamages


*Count the number of zip codes with solar installations under the reallocation scenarios
count if !missing(num_systems) /*original*/
count if state_new_num_systems>0
count if state_new_num_systems30>0

*Calculate mean avoided damages across zip codes with installations per avg system
mean total_damages if !missing(num_systems) /*original*/
mean total_damages if state_new_num_systems>0
mean total_damages if state_new_num_systems30>0

*Calculate the avoided damages (and the change in avoided damages) by state
bysort state (rankbystate): egen sumTDbyState = sum(newState_total_damages)
bysort state (rankbystate): egen sumTDbyState30 = sum(newState_total_damages30)
bysort state (rankbystate): egen sumTDbyState_Original = sum(total_damages_all)
gen state_gain = (sumTDbyState - sumTDbyState_Original) / sumTDbyState_Original
gen state_gain30 = (sumTDbyState30 - sumTDbyState_Original) / sumTDbyState_Original

save `"${data_intermediate}/Simulation_Results_Intrastate"', replace
outsheet zipcode NAME state POPULAT POP_SQM pcac INT total_damages num_systems state_new_num_systems state_new_num_systems30 using `"${data_intermediate}/Simulation_Results_Intrastate.csv"', comma replace

collapse sumTDbyState sumTDbyState30 sumTDbyState_Original state_gain state_gain30, by(state)
rename sumTDbyState damages_new
rename sumTDbyState30 damages_new30
rename sumTDbyState_Original damages_old
outsheet state damages_new damages_new30 damages_old state_gain state_gain30 using `"${data_intermediate}/state_reallocation_avoided_damages.csv"', comma replace

save `"${data_intermediate}/state_reallocation_avoided_damages"', replace

